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Abstract 

We analyze the phase equihbria of systems of polydisperse hydrocar- 
bons by means of the recently introduced moment method. Hydrocarbons 
are modelled with the Soave-Redlick-Kwong and Peng-Robinson equations 
of states. Numerical results show no particular qualitative difference be- 
tween the two equations of states. Furthermore, in general the moment 
method proves to be an excellent method for solving phase equilibria of 
polydisperse systems, showing excellent agreement with previous results 
and allowing a great improvement in generality of the numerical scheme 
and speed of computation. 

1 Introduction 

In this paper we analyze the phase behaviour of a mixture of hydrocarbons, 
by means of the moment method [3 |H]. This method allows to reduce the 
number of degrees of freedom of the free energy, which normally depends on the 
concentration of each specie in the mixture, to a smaller number of moments 
of the density distribution which already appear in the excess part of the free 
energy. By doing this, one is able is reduce the number of equations needed 
to analyze the phase equilibria and, at the same time, by projecting the free 
energy onto the space generated by the moments only, to check for global and 
local stability of the phases [5]. 

The approximation made when introducing the moment free energy can be 
efficiently controlled and minimized by means of the adaptive method of choice 
of extra moments which allows to reduce the deviation of the moment 
method solution from the exact solution, by simply retaining two extra moments, 
beyond the ones appearing in the excess free energy. This iterative method, 
which it can be proven, converges to the exact solution, as long as it converges at 
all, shows to represent an excellent compromise between approximation, which 
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can easily be reduced to an error smaller than 0.01%, and computational speed. 
Furthermore, the resulting algorithm turns out to be stable and to be very 
little affected by the number of species in the mixture. In fact, as the number 
of unknowns is not increased by the increase of the number of species, the 
computation is hardly affected at all, with just a small influence on its global 
speed of computation, while no relevant error is introduced. 

The numerical results agree very well with the results obtained with a widely 
diffused commercial program in different points of the phase diagram. The 
concentration of each component in the coexisting phases and the density of both 
phases are evaluated correctly. Clud point is detected exactly. Furthermore, the 
introduction of heavy species, up to n-C15, even present in very small amount, 
does not compromise either numerical results or computation. 

2 Polydisperse hydrocarbons 

In order to analyze the phase equilibria of hydrocarbons, we will refer to the two 
equations of state most widely used to describe them, i.e., the Soave-Redlick- 
Kwong (SRK) |6j equation of state and the Peng-Robinson (PR) [5] equation of 
state. Both are cubic equations of state and thus are able to predict gas-liquid 
phase transitions. Although originally introduced for pure systems, as we will 
see, they are both easily extended to describe multicomponent,i.e., polydisperse 
systems. As we will show, given the polydisperse form of two equations of 
state, one can easily obtain the Gibbs and Helmoltz free energies, by Legendre 
transforming, and therefore obtain the phase equilibrium equations that are to 
be solved, in order to fully analyze the phase behaviour of the system. 
The SRK equation of state is generally written as 



where N is the total number of particles, V the total volume, kb the Boltzman 
constant and the parameters a, b and a(T) depend on the critical temperature 
and pressure, shape and size of the molucules etc., of the specific hydrocarbon. 

The extension of the equation above to the case of polydisperse system is 
rather straightforward, if one introduces a set of parameters ac^i, hi, ai(T) for 
each specie i and defines new global parameters B and D as follows 



P = 



NkbT N^a{T)ac 
V ~b ^ V{V + b) 



(1) 




(2) 



and 




(3) 



where 




(4) 
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In this way, the polydisperse version of Eq. ([T]) is 

p N.^T D 

V~B V{V + B) ^ ' 

where = Ni is the total number of particles of the system. 

The Helmoltz free energy can now be obtained simply by solving the equation 

dV 

and by introducing F\a-, the ideal part of the free energy, i.e., the free energy of 
a mixture of ideal gas 

i^id = 'tBT fin ^ - l) (6) 



The free energy then turns out to be 



V D V + B 
F(n, y, T) = Fid + NkbT In ^-^ - - In (7) 

The above quantity is extensive, one can therefore define an intensive "free 
energy density" / = F/Vol. Introducing a density distribution p{k) = N^/V 
and multiplying by /? = I/k^T, the free energy density turns out to be 

/3/[p(fc),T] =/3/id-pln(l-B)-£ln(l + B) (8) 

B 

where the ideal part is just 

/Ki = A^BT^p(fc)[lnp(fc)-l] (9) 



and we have defined two new parameters B and D, by rescaling B and D with 
the volume 

^ = ^ = J2bkPik) (10) 



V 

^ = ^=^/3a,,,(r)p(fc)pO-) (11) 



The non-ideal part of the free energy in Eq. © is called excess free energy / and 
contains the terms due to the interaction between the particles in a non-ideal 
gas. 

The Gibbs free energy can now be obtained from the expressions above, 
simply by Legendre transforming F as G{N, P, T) = miny{F(A^, V, T) + PV}. 
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We now introduce a Gibbs free energy per particle g = G/N and divide again 
the resulting function into the ideal part 

/35id[x(fc),P,T] = ^ = ^x(fc)lna;(fc) +ln/3P (12) 



and the excess part 



/35[x(fc), P, T] = ^ = - V x{k) [In + 1] + ^ + ^ (13) 

k P P 

where the number fraction x{k) of the specie k is just x{k) = N^/N = p{k)/ p, 
with p = p{k) = A^/y the overall density. 

From the above equations (jl2ll3p wc can now derive the phase equilibrium 
equations pf. = p\, for the coexisting phases a and 6, and each specie k as 
Pk = dG/dNk — dg/dx{k). For a system of M species, dividing in P phases, 
the phase equilibrium is therefore fully analyzed by solving a system of (P— 1)M 
equations, plus the M equations given by the conservation of the total number 
of particles, i.e., X]a^"(^) ~ x'^^Hk), where x^^\k) is the number density of the 
fcth specie of the parent, in the MP unknowns x°'{k). The values of P and T 
arc set as external parameters. 

As far as the Peng-Robinson equation of state is concerned, not much change 
is needed in the equations above. The PR equation is generally written as 

P=^ ^ (14) 

V-B V{V + B) + B{V - B) ^ ' 

where B and D differ from the SRK case in the numerical coefficients of and 
ac,k- Once again, from the equation above one gets the excess part of Helmoltz 
free energy as 



/3/ = P In + V 1^ In TT^-r^^ (15) 




Similarly, the excess part of the Gibbs free energy per particle turns out to be 
/35 = -^.T(fc)[ln^PF + l] + ^ + ^ (16) 

2.1 Truncatable systems 

A polydispersc system is said to be truncatable when the excess part of its free 
energy, say the Helmoltz free energy, is a function a limited number of moments 
Pi of the density distribution /^(fc) 

p,=Y.w,(k)p{k) (17) 

k 
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with given weight functions Wi{k). In other words, the Hehnoltz free energy of 
a truncatable system is 



/[p(fc),T] = ^ p{k) [Inp(fc) - 1] + fip,) 

k 

For a truncatable system one has 

= §^ - Inp(fc) + Mk)PfL^ 

where the excess moment chemical potentials pn = df / pi. By imposing the 
equality of the chemical potentials in all the coexisting phases, one gets that 
the density distribution of the coexisting phases must have the form 

p'^(fc) = R{k)eM-PY.~^>^m (18) 

i 

By imposing the lever rule, i.e., the conservation of the total number of particles 
per specie, '^a''^"' P°'(^) — P^^H^) is the density of the parent and v°- = 

V^/V is the volume occupied by the phase) one finds that the function R{k) 
has the form 

mk) . = (19) 

^ ' Ea«"exp[-/?A-(fc)] Ea«"CXpME./i>.(fc)] 

Although formally solved through the two equations above, an actual numerical 
solution of the system is not easily found. Eq. actually represents a set of, 
say, M (for M species of particles) self consistent all strongly coupled through 
the denominator of Eq. ([T9| . The problem is that, although the excess free 
energy is a function just of the moments pi, usually a much smaller number 
than the number of species, the ideal part of the free energy is still function 
of the whole density distribution p{k). Ideally, one would like to reduce the 
problem to a smaller number of degrees of freedom, by expressing also the ideal 
part of the free energy as a function of the moments only. While this argument 
will be treated in the next section, here we will show that both the SRK and 
the PR equations of state are in fact truncatable. 

In fact, it is rather easy to show that the equations of state generate two 
truncatable systems if one introduces two moments of the density distribution 
pi and p2 as follows 

p, = B = Y,bkP{k) (20) 

fe 

P2 = Y.'^^P^^) (21) 



where, from Eq. (|4]), = ^ /3ac\kCtk{T). From the definition above and Eq. ([3]), 
we get that D = p\. Plugging (|20l2ip into Eq. ([S]), we therefore get, for the 
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SRK equation of states 

Php.P1.P2) = -pln(l - pi) - ^111(1 + (22) 

Pi 

Note that in the above expression, the overaU density p is itself a moment of 
the density distribution, with weight function woik) = 1, as p = = X^A; Pi^)- 
In other words, the excess part of the free energy is fuhy described by the 
knowledge of only three moments of the density distribution po,pi,p2 and not 
on the whole distribution p{k) itself. 

For the PR equation of state the argument is again similar to the case of the 
SRK equation of states. Introducing again the three moments pq, pi and p2, 
we get that the excess part of the Helmoltz free energy is just 

/3/(po, pi,P2) = -p\n{l - Pi) - -. In , , (23) 

4 pi \^1 + (1 + V2)pi J 

As far as the Gibbs free energy is concerned, it is easy to show [8] that the Gibbs 
free energy inherits its moment structure from the Helmoltz free energy. How- 
ever, this time, one usually introduces normalized moments rrii of the number 
density distribution x{k) defined as 

m, = y w,(fc)a;(fc) = (24) 

Clearly this time mo = X]fe^(^) ^ thus, since / depends on three moments, 
Po, Pi, P2, the Gibbs free energy depends itself on the overall density, which, 
however, is obtained from the equation of state, as a function of P, mi and m2. 
In other words, the Gibbs free energy turns out to have one degree of freedom 
less than the Helmoltz free energy. 

In any case, with the definitions above, one gets that the excess part of the 
Gibbs free energy for the SRK equation of states is just 

/3g(mi,m2) =ln|^ -1 + /?— -ln(l-pomi)- ^ln(l+poTOi) (25) 
pF po nil 

Similarly, for the PR equation of state, we get 

R~f \ 1 PO wi A^V2m2 I' i + (i- V2)pomi ^ 

^.g(mi,m2) = In — + lnl-pomi) + - In y=- 

PP Po 4 mi \^ 1 + (1 + V2)pom2 / 

(26) 

Note that in the two equations above we have omitted the dependence of g on 
P and T. 



3 The moment method 

Truncatable systems allow to express the excess free energy as a function of a 
small, say M, number of moments only. However, as we saw in the previous 
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section, the difficulty of solving the phase coexistence equations remains largely 
unaltered, as the ideal part of the free energy is still function of the full density 
(or number) distribution. Ideally, one would like to express the ideal free energy 
too, as a function of the moments only. This is in fact possible, by means of the 
moment method [H [71 [S] . While the following description will refer mostly to 
the Hclmoltz free energy, similar considerations can be made for the Gibbs free 
energy [8], with the introduction of the normalized moments and the number 
density distribution. 

The moment method arises from the hypothesis, in fact verified in different 
works [Tl[71[TTJ[Tn|, that the excess free energy is mostly responsible for the phase 
behaviour of the whole system. This is in fact not surprising, as the ideal free 
energy is overall convex, and thus does not allow for phase separation. With 
this in mind, the moment free energy is constructed as follows. We subtract 
from the actual free energy a term p{k) In p'^^^k) , where p'^°-'(fc) is the density 
distribution of the parent. This term, as linear in the density p{k), docs not 
affect the phase behaviour, as it adds just a constant to the chemical potential 
fj.{k) = df /dp{k). The resulting function is then minimized with respect to p{k) 
with the M moments appearing in the excess part as constraints (2 in the two 
previous cases). The minimum value of the resulting free energy is then found 
to be 

M 

i 

where p is just a vector having the moments pi as components and the As are 
the M Lagrange multipliers. The minimum value of the free energy is reached 
for a density distribution from the family 

Pmom(fc) - P^°^ (fc) exp ^ X^w^ik)^ (28) 

From the moment free energy (|27p . one can also define the moment chemical 
potentials /i^, as fii = dfnwm/dpi = \i+ jxi. The pressure is obtained from the 
Gibbs-Duhem relation as 

p = Y^p{k)p{k) ^ f ^Y.^'^p^- f 

k i 

It is easy to show that the above expression obtained from the moment free 
energy is in fact identical to the one obtained from the exact free energy [5]. 
Furthermore, it is easy to show that the moment free energy correctly detects 
the onset of phase coexistence. In other words, one can show [51 [S] that any 
two phases coexist for the full system, if and only if, they coexist for the mo- 
ment free energy, thus, p.^{k) ~ p''{k) <^ p,°; ~ p'^. Thus, at least up to the 
onset of the phase coexistence, the full solution in Eq. ([T5)) actually belongs to 
the family in Eq. ([28)) . Cloud point and shadow phases are then correctly de- 
tected by the moment free energy solution above. Furthermore, spinodals and 
critical/tricricital points arc found exactly [S]. 
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The enforcement of the lever rule only for the moments is in fact the only 
approximation we make in using the moment method, as this does not ensure 
the satisfaction of the complete levere rule, while, as mentioned earlier, equality 
of pressure and chemical potentials are ensured. However, as shown in details 
in [TJ 1111 110] , the approximation can be reduced efficiently by retaining extra 
moments and, in particular, by means of the adaptive method of choice of extra 
weight functions which allows to obtain a solution as close as wanted to the 
exact one, by retaining only two extra moments. 

In order to give a more precise insight of the actual problem one has to solve 
in the case of the two equations of state mentioned earlier, let us sketch the 
resulting system of equations, obtained within the moment method approach. 
As mentioned, since we have P and T as external parameters, we move on to 
the Gibbs formalism. Thus we evaluate the gibbs free energy and from that, 
we calculate the moment chemical potentials as = dgmom/drrii, where rrii are 
the normalized moments. Let us now assume we have a Gas-Liquid demixing 
and let us call = N°'/N^°'> the fraction of particles in each phase (G or 
L). The phase coexistence is then fully solved by enforcing the equality of the 
moment chemical potentials and of the quantity 11 = gmom — In P — X^i^^o "^iMi j 
which is a sort of Legendre transform of the pressure [5], in all the coexisting 
phases. We must also enforce the conservation of the total number of particles, 
i.e., J2a^'^i^) = N'^^^k). If we multiply by Wi{k) on both sides and sum 
over k, we get, after rearranging, the lever rule for the normalized moments 

(f)°-m^ = rnf \ which is the condition we actually enforce. Thus, for the two 
equations of state, the system of equations we have to solve turns out to be 

mf^ = (fPm^ + (1 - (l)^)m\ 

m^2^ = (fPm'^ + (1 - 0'^)"i2 
P = P(p^,mf,mG) 
P = P(pL,mL,mL) 

i.e., 7 equations in the 7 unknowns Xf , X\ , , , Pqj 4>'^ ■ The moment 

chemical potentials and the pressure are then for the SRK equation of state: 

^^1 = PXi + - — -H ^\n{l + pomi) 

I — pomi ■mi[l + pQini) ml 

/3/i2 = /?A2-2^ ln(l + pomi) 

771 1 



(3P = 



Po Plml 



1 - pomi 1 + poTTii 
For the PR equation of state, the calculations are just slightly more complicated. 
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as now we get: 



/3P 



/3A2 



1 - po?TM 
V2 m2 



2 mi 



In 



mi(l + 2pomi - plml) 
1 + (1 - V2)pomi 



V2 ml , 
— 7T In 



1 + (1 - V2)po?71l 

1 + (1 + V2)pQmi 



1 + (1 + V2) 



plml 



4 Numerical results 

As a way of example, in order to show the potential of the method described 
in the previous section, in this section we will show some numerical results 
obtained by applying it to a real case-study. We solve the phase equilibria for 
a mixture of 24 hydrocarbons, up to nC15, along a straight line crossing the 
{P, T) phase diagram. The calculations are done using the SRK equation of 
state, although no relevant difference is observed when using the PR equation 
of state. Our numerical results are compared with the results obtained with a 
commercial program (PVTsim), licensed to Snamprogetti s.p.a. (that provided 
the results). In Fig. SI we plot the concentration in mole % of CI (methane), 
in the two coexisting phases, against the temperature. Our result agrees very 
well with the points obtained with PVTsim. The concentration of both phases 
is evaluated correctly. Furthermore, the cloud point, i.e., the point at which 
the liquid phase appears, is detected exactly on the phase envelope shown by 
PVTsim. 

In Fig.|4]we show again the same case. Now we plot the concentration of n-C4 
against the pressure. As the pressure-temperature path enters the coexistence 
region, liquid is found. This time the component represents less than 1% of the 
total composition of the gas, while it is about 10% of the total composition of 
the liquid. Again, even with a heavier hydrocarbon, our numerical results are 
in excellent agreement with the ones obtained with PVTsim. 



5 Conclusion 



We have applied the moment method to the analysis of phase equilibria of mix- 
ture of hydrocarbons, using the SRK and PR equations of state. Our results 
show that not only the moment method is applicable as it correctly detects 
gas-liquid phase coexistence. Even with a large number of components in the 
mixture, our algorithm remains robust and numerical calculation fast. Fur- 
thermore, our numerical results agree quantitatively very well with the results 
obtained using a widely used commercial program (PVTsim) . No further demix- 
ing, beyond the G-L coexistence is observed, however, it is not yet clear whether 
this depends on the equations of state used, or rather to the choice of the mix- 
ture of hydrocarbons. It may be possible to have the coexistence of more than 
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Figure 1: Concentration of methane against temperature in the two coexisting 
phasese (sohd and dashed hues), compared to the resuhs obtained with PVTsim 
(diamonds). The hquid phase appears correctly at lower temperature as the 
path across the phase diagram crosses the phase envelope. The cloud point is 
detected exactly by our method. Deep inside the coexistence region some small 
deviations appear, although it is not clear whether they are due to the moment 
method, or to the approximations introduced by PVTsim. 

two (G and L) phases, e.g., more than one different liquid and/or gas phases, 
using a wider distribution of hydrocarbons. 

As far as future developments are concerned, one could proceed to the anal- 
ysis of a fully polydisperse case, i.e., by introducing a continuous distribution 
of species. Clearly the two equations should first be extended to the continuous 
case, e.g., by introducing a continuous dependence of the acentric factor on the 
size of particles. It is possible that further demixing appears using different 
distributions. The extension to the continuous case should not have any impact 
on our method. Clearly, the computation may be slightly slower, as in that 
case the moments have to be evaluated by integration rather than summation 
over a finite number of species; however, no formal or substantial adjustment is 
needed. 

Finally, further applications of the moment method could be tried. For 
instance, one could try to apply it to metallurgy or other cases of industrial 
interest. Clearly, when solid phases appear, one would have to extend the 
analysis introducing spatial and/or orientational degrees of freedom. However, 
as previous results already show [TJ [31 [31 SI [TI] [TU] this should not represent a 
limitation for the method. 



10 



mphase G 

PVT G + 

mphase L — 
PVT L X 




Figure 2: Concentration in mole % of n-C4 against pressure across the phase 
diagram. As pressure and temperature drop enough to enter the phase coexis- 
tence region, the liquid phase is found. The concentration the element in both 
phases are in excellent agreement with the ones obtained with PVTsim. Again, 
it is not clear whether the small deviations inside the coexistence region are 
due to our method, or rather to approximations and truncations introduced by 
PVTsim. 
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